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Abstract 

We present a finite element technique for the efficient generation of lower and upper 
bounds to outputs which are linear functionals of the solutions to the incompressible 
Stokes equations in two space dimensions: the finite element discretization is effected by 
Crouzeix-Raviart elements, the discontinuous pressure approximation of which is cen- 
tral to our approach. The bounds are based upon the construction of an augmented 
Lagrangian: the objective is a quadratic “energy" 1 reformulation of the desired output: 
the constraints are the finite element equilibrium equations (including the incompress- 
ibility constraint), and the intersubdomain continuity conditions on velocity. Appeal to 
the dual max-min problem for appropriately chosen candidate Lagrange multipliers then 
yields inexpensive bounds for the output associated with a fine - mesh discretization, the 
Lagrange multipliers are generated by exploiting an associated coarse-mesh approxima- 
tion. In addition to the requisite coarse-mesh calculations, the bound technique requires 
solution only of local subdomain Stokes problems on the fine-mesh. The method is il- 
lustrated for the Stokes equations, in which the outputs of interest are the flowrate past, 
and the lift force on, a body immersed in a channel. 


1 Introduction 

Fast solvers are essential in engineering design due to the large number of appeals to the 
simulation performed within a design cycle. Indeed, the search for faster solution strate- 
gies remains a major research objective. Parallel computing, domain decomposition, pre- 
conditioners, higher-order schemes, and adaptive methods are just some of the successful 
techniques that are being brought to bear on this important problem. 
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In actual practice, a typical design effort consists of the optimization of an objective 
function with respect to selected design variables. The quantity of interest in such objective 
functions is typically not the entire field solution, but rather a characteristic metric of the 
system that we will term an “output”: for design applications, this output value is more 
relevant than the entire field solution. Recently, a fast approach has been developed [19, 20, 
21, 22] to calculate rigorous bounds to an output at a fraction of the cost of a traditional 
solver. 

This technique calculates lower and upper bounds to outputs which are linear functionals 
of the solution to coercive partial differential equations; recent extension to noncoercive and 
nonlinear problems is discussed in [25]. The bounds are for the output associated with a 
verv accurate spatial discretization which we shall call the truth mesh, direct calculation 
of the output on this discretization would be extremely expensive. In our approach, the 
computation of the bounds nevertheless remains inexpensive; consisting of only global solves 
on a coarse— mesh: domain decomposition performed along the edges of this coarse— mesh, 
and finallv, calculations of local subdomain Neumann problems on the truth mesh. In fact 
the coarse-mesh may be considered as the “working” mesh utilized in a design cycle: the 
bound values then serve to relate the accuracy of the design optimization to the “truth". The 
technique is based on the construction of an augmented Lagrangian, in which the objective 
is a quadratic energy reformulation of the desired output, and the constraints are the finite 
element equilibrium conditions and inter-subdomain continuity requirements; the bounds 
are then derived by evoking the dual max-min problem for appropriately chosen candidate 
Lagrange multipliers. 

In this paper we extend this technique to the incompressible Stokes problem [21], of inter- 
est in its own right, but also as a precursor to the incompressible Navier-Stokes equations. 
The new considerations addressed are threefold. First, the Stokes problem is itself a con- 
strained minimization problem. Therefore our Lagrangian must be modified to include an 
additional primal variable, the pressure, and an additional Lagrange multiplier to impose 
the incompressibility constraint. Second, the pressure term contained in our new Lagrangian 
will not be controlled by the energy term, which may thus lead to infinite bounds. The so- 
lution to this difficulty is the use of the Crouzeix-Raviart element [15, 26], which permits 
us to exactly eliminate the dependence of the Lagrangian on the pressure variable through 
a projection technique which, thanks to the discontinuous (and hence decoupled) pressure 
space, can be effected solely through problems local to each element. Third, higher order 
velocity approximation is required in the Stokes problem to satisfy the inf-sup condition. 
This requirement also necessitates higher-order hybrid flux construction, which is developed 
here in a formulation similar to that described in [5]. Finally, regarding computational sav- 
ings, the domain decomposition of the Stokes problem offers even more substantial savings 
than the domain decomposition of elliptic problems. 

Our work has benefited from previous efforts in the a posteriori error estimation commu- 
nity [3, 5, 7, S, 9, 10, 13, 14, 17, 27, 30]. In earlier papers [19, 20, 22] we have described 
the similarities between our bounds technique and both “explicit” and “implicit” a poste- 
riori error estimators for elliptic equations. Similarly to earlier implicit techniques, we base 
our bounds technique on local independent subproblem calculations. However, our bounds 
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offer the advantage of measuring the error in norms different than the energy norm. Indeed, 
for quantitative confirmation of engineering design quantities, we must measure the error 
directly in the norm associated with these outputs. Recent explicit error indicators for the 
error in linear-functional outputs have been developed [14] based on the Aubin-Nitsche du- 
ality procedure. These techniques allow adaptive improvement of finite element predictions 
for the desired engineering output, and can also be applied to the Navier-Stokes equations. 
However, these estimates are less quantitative than ours due to the presence of constants that 
cannot be precisely evaluated, and thus the goal of design confirmation is less satisfactory 
achieved. 

For the Stokes problem, Verfurth [31] has developed implicit error estimates based on the 
solution of local Stokes problems, and explicit estimators based on a suitable evaluation of 
the residual of the finite element solution, that provide estimates for the error in the energy 
norm. In [12], Bank and Welfred successfully reconsider the implicit error estimators for 
the Stokes problem. A comparison of all of these methods [11] indicates that the estimates 
are a good indicator of the error, that the explicit estimator is about a factor of two less 
expensive than the implicit estimators, and that the implicit estimators require about one 
fourth the computing time needed for the solution process. We note that these estimators 
have been developed for the mini-element discretization of the Stokes problem [2], which is 
based on piecewise continuous linear velocities augmented with quadratic bubble functions 
and piecewise continuous linear pressures. Our technique is limited to the discretization of 
the Stokes problem by Crouzeix-Raviart (discontinuous pressure) elements. 

Less standard approaches to measuring the error have been proposed in [18. 4]. The 
error estimators proposed by Ladeveze et al. [18] measure the error in the constitutive law 
of materials in the limit of incompressible solids. (Recall that there is a direct analogy 
between an incompressible linear-elastic isotropic solid in equilibrium and an incompressible 
Newtonian fluid in the steady creeping limit [23].) Another implicit estimator for the Stokes 
problem is found in [4]: in this approach the error estimator is based on local residual 
problems which require only the solution of decoupled subdomain problems of Poisson type 
with Neumann data. Although this method has the advantage of being faster than implicit 
methods that require solution of local Stokes problems, the bounds obtained are for an 
“equivalent 7 ’ energy norm, and thus not directly relevant to validation and confirmation in 
engineering design. 

We remark that most of the previous work on a posteriori Stokes error analysis is focused 
on estimating the error for application to mesh adaptivity rather than directly addressing 
engineering design problems; there is a relative lack of methods for validation and confirma- 
tion which focus on rigorously quantifying the error in the outputs of interest. Nevertheless, 
the utility of adaptive mesh technology indicates that our technique must be extended to 
quantify the error locally for use in adaptive error control procedures. For elliptic Partial 
Differential Equations such an extension has already been presented in [24]; generalization 
to the Stokes problem, though not considered here, should be relatively straightforward. 

The outline of the remainder of the paper is as follows. In Section 2 we describe our model 
problem and the output linear functionals which we will investigate. In Section 3 we present 
the finite element discretization to which we apply our bounds technique. In Sections 4 and 5 
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we describe and prove the bounds procedure and the properties of our estimates. In Section 
6 we develop an approach to decrease the bound gap. Finally, in Section 7, we illustrate 
our technique for the Stokes problem and associated outputs of interest to demonstrate the 
engineering relevance of our technique. 


2 Model Problem 

2.1 Governing Equations 

We consider the steady creeping flow of an incompressible (p = constant) Newtonian fluid 
with constant dynamic viscosity, p, between two plates with a periodic array of rectangle 
obstacles in the center. This geometry is presented in Figure 1. where (x 1 ,o- 2 ) denotes the 
coordinate system with corresponding unit vectors ej, e 2 ; is the domain; and Tj, j = 
are the domain boundary segments. The flow is driven by a forcing term which can 
be interpreted as a pressure gradient A P/L in the ei direction. The velocity and pressure 
perturbations are periodic in the ei direction. 

To describe this flow we use the “Laplacian" form of the incompressible Stokes equations. 
In Gibbsian notation, the velocity vector u and the scalar perturbation pressure field p satisfy 

— Au + Vp = f, in H, (1) 

— V • u = 0, in (2) 

with no-slip Dirichlet and periodic boundary conditions, 

u = 0 on T;, i € {1, 3, 5}. (3) 

u lr 2 = u lr 4 , ( 4 ) 

Vu|r 2 = Vu|r 4 . (5) 

Here f is the volumetric force, f T = [10]; for convenience we set the viscosity to unity. We 
also require that Jq p d.4 = 0 for uniqueness. 

The variational form of (l)-(2) is: Given f € (7"t _1 (ff)) 2 , find u = (ui,u 2 ) € CHj(n)) 2 


and p <E Lq)!)) suc h that 

[ Vv- Vu-pVv-vf dA = 0 VvGHj(fl)®Wj(fl), (6) 

Jn 

-fqVudA = 0 VqeLl(Q), (7) 

Jq 

where d^4 is a differential area element, and 

=v|r 4 ; v|r, = 0, i€ {1,3,5}}, (8) 

I 2 (H) = {qeL 2 m [ qdA = 0}, (9) 


where Ti 1 ^) and L 2 (Q) are the usual Sobolev spaces [1]. We also introduce X = 7i l (Vl) ® 

WHH), y = I 2 (fi), x° = wj(ft) ® K(Q), and y° = L 2 (Q). 
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Figure 1: Computational domain: Ti, T 3 and r 5 are homogeneous Dirichlet boundaries, 
and r 2 are periodic boundaries. 


2.2 Output Linear Functionals 

We assume that our output, s, may be expressed as a linear (or more generally, affine) 
functional of the velocity components u. and a linear functional of the pressure p, that is 
s = C(u ,p) - £ v ’(u) + £ F (p) where 


C : X © y -+ R. 


( 10 ) 


or 


£ v : X — ► R, (11) 

t p R. 


It is clear that t is a linear functional on the product space X x y. On physical grounds, 
£ p (l) = 0, since the pressure level is arbitrary, and thus must not affect the output; the 
mathematical ramifications of this condition will become clear later. 

Examples of possible linear functionals include the flowrate through the channel, or the 
lift force on the body immersed in the fluid. The particular linear functional for the flowrate 
(output s* 1 ') is defined as 


£ v '(v) = ^ / v • ei dA, Vv € X, 
U Jn 

t p (q) = 0, Vg € >\ 


( 12 ) 


where L(= 2) is the height between the plates. Note that this output functionals is bounded 
for all v in X. Another important engineering output of interest (s (2) ) is the lift force acting 
on a body. We evaluate this force with the following functionals: 

^ V ( v ) = / Vx-Vv-x-fd.4, 

Jn 


5 




e F M = - I qv-x < 1 - 4 , 

Jq 


( 13 ) 


or equivalently 


„( 2 ) 


= / V X 
in 


Vu — p'V ■ x ~ X ' f dA 


(14) 


where x is an y continuous function in X such that x ' e 2 — 1 on Ts and x — ® on ^e other 
non-periodic boundaries. 

The motivation behind the choice (13) is once again to obtain bounded functionals, since 
we can easily predict specific convergence properties only for £ v € ( 0 ) and i p € L 2 (Cl)\ 

it is shown in [21] that (13) is indeed bounded. To show that we correctly reproduce the lift, 
we first note that it corresponds to 


5 ( 2 ) _ / V • (x • Vu) + V • ((x • V)u) - V • (p X ) d.4 - / V • ((x ■ V)u) d.4. (15) 

in 

By application of Gauss’ theorem we then obtain 

s (2) =/ x ' (<7n)ds — / V • ((x • V)u) d.4, (16) 

Jc> n in 

where a is the stress tensor and n is the outward normal vector on the domain boundary. 
Finally, we demonstrate in [21] that the term f n V-((x-V)u) d,4 is zero for smooth solutions, 
since both the tangential and normal derivatives of the normal velocity vanish, the latter 
thanks to incompressibility; from our definition of (16), s (2] thus reduces to the lift. (Note 
that future work should consider the stress formulation of the Stokes equations, which more 
naturally generates the stress contributions on the boundary.) The functional. (13), also 
permits the calculation of the drag force acting on the body similarly bv choosing a function 
X such that x ' e i = 1 on IV 

We close this section with two remarks. First, if we choose x to be an incompressible 
field, then the pressure part of the functional {i p {p)) is zero. To show that this choice is 
compatible with the boundary condition, we apply Gauss' theorem to find 

f V • X = f X ■ n ds = f x • n ds = 0, (17) 

Jq JdQ J r 5 

where the final equality obtains since T 5 is a closed boundary contour and x • e 2 = 1 on r 5 . 
Second, we note that (17) also proves that, as required, ^ p (l) = 0 in (13). 


3 Finite Element Discretization 

We first introduce the necessary triangulations, and the general finite element ingredients, 
such as the bilinear and linear forms and function spaces, that will be required in what 
follows. 
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3.1 Triangulations 

Two different types of triangulations are required for our “hierarchical” bound procedure, 
the H - mesh and the h- mesh, where the latter is a refinement of the former. The h - mesh 
is our fine mesh, which serves as the “truth” mesh; by “truth” we indicate our assumption 
that the difference between the numerical solution obtained for this fine-mesh and the exact 
solution is negligible. The H - mesh is our “working” mesh, which is used, in conjunction 
with local Stokes problems, to calculate the bounds. 

As our H - mesh discretization of Q we take a geometrically conforming regular triangu- 
lation Th consisting of K triangles Th such that 

n= u t h . (is) 

ThZTh 

We denote the set of all (open) edges 7 of this triangulation as £(T//), and the set of three 
edges 7 associated with each element Th as £(Th)- We denote the set of interior edges 
as S int (T„), and the sets of Dirichlet edges — the edges that, are part of Dirichlet boundary 
segments — as £ q ( Th ) - We denote the set of A vertices of the triangulation by *\A(Tfj). 

The triangulation and elemental edges are, of course, related. In particular, given an 
edge ~jj H in £{Th)- we shall indicate the coincident edge 7 in £{Th) as 7 = E(~/t h )- We 
next associate with each edge 7 in £(T H ) a unique normal n 7 such that, if 7 lies on <9D. n' 
coincides with the outward normal n on <9fi. Then, for all Th in Th . and all edges 7 r H in 
£(Th), we define 

aj T H H = h EhT » ] -n> T n, (19) 

where n 7T H is the outward normal on 7T„ with respect to T H ■ In essence, a?* is ±1 on the 
two “sides” of an edge 7 in £(Th)- 

Finally, we introduce the h - mesh triangulation T h consisting of triangles T h such that 

n= |J r h . (20) 

Th€'Th 

We shall require that Tj bea refinement of T H , in that we can express each T H in T H as 

T h = U T h , (21) 

where 7 Z Th is thus the set of h - mesh elements that comprise T H . A uniform R refinement 
will denote an h - mesh in which 7 Zj H consists of R 2 triangles Th similar to Th- 

3.2 Bilinear and Linear Forms 

We define here the bilinear and linear forms required for the Stokes problem. We first need 
to define a “broken” space in which no continuity is required on between elements; this space 
serves to define functions on the local subdomains. In particular, we define 

«l(ft) = {ve L 2 m v\ Th G H\Th), VTh G Th }, (22) 
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and associated product spaces X' — ('Hj(fi)) 2 . Xr H — ('H 1 (Th)) 2 , and Tt w — L 2 (Th)- 
We now define the bilinear form associated with our operators as 



a(w,v)= a T H (w\T H ,v\T H ), V(uj.t’) € (Wi(n)) 2 , 

Th^Th 

(23) 

where for all Th in Th 



a TH (iu, v) = [ Vw Vv d.4, V(tr, v) G (H x (T h )) 2 . 
Jt h 

(24) 

In addition 

we denote 



a(w,v) = a(u7j,t>x) + a(iu 2 ,v 2 ), V(w.v) G (-T*) 2 . 

(25) 

and 

ar H { w.v) = ar H (iUi,ri) + c, t h ( U' 2 ■> v 2 ) ■ V(w.v) G {X Th ) 2 ■ 

(26) 

Similarly, 

d(w,q) = ^t h ( w \t h - q\T H )- V(ic, q) G 'Hl(H) G L‘(Q.). 

Th6T h 

(27) 

where for all Th in Th 



di r (w,q)=[ qVweidA . V(uu g) € ^(T w ) ® i 2 (ft), 

Jt h 

(28) 

and 

dr H (w, 9 ) = 4»( w i»9) + 4 h ( U! 2’ 9), v (w, 9) e ® 3V* • 

(29) 


Note that a and d correspond to the Laplacian and divergence operators, respectively. 

We next introduce a set of "jump'’ bilinear and linear forms required in our variational 
formulation. These forms will be applied in a scalar fashion to each component of the 
velocity. In particular, we define the bilinear from 

b(w,t)= £ £ J w\ Th t\ EhTH) ds, VOM) € «I(n) X e, (30) 

and 

6(w,t) = 6(117, ti) + b( W 2 ■ <2), V(w,t) G X‘ x Q 2 , (31) 

where iv\ j H in (30) is to be interpreted as the trace of w\j H on 7 r H . and Q = H~ x ^ 2 {S{Th))\ 
note that t is defined only over the edges of the triangulation. Effectively, (30) computes the 
moments of the jumps in w over internal edges, and the moments of w over boundary edges. 

We now introduce our linear functionals. Associated with the volumetric inhomogeneities 
we have 

i N {w) = H Vw € X\ (32) 

where for all Th in Th 

£j h (w) = f w • f di4, Vw € Xt h - (33) 
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Associated with our output functional, we introduce 

<°'(w) = E d( w lr„), Vw 6 (34) 

Th^Th 

such that 

^°"(w) = ^(w), VwGA\ (35) 

Similarly we can introduce a linear functional for the pressure, 

i OP (q)= E P’^Th): VgGjF 0 . (36) 

Th^Th 

such that 

e 0P {q ) = * P (?), V 9 € y\ (37) 

Here l v () and i p {) are the formal output functionals introduced in ( 11 ). 


3.3 Function Spaces 

As already indicated, we consider two different spatial discretizations: 8 = H and 8 = 
h , which correspond to our “working” and “truth" discretizations, respectively. For the 
Crouzeix-Raviart approximation spaces of interest [16, 26] the velocity space is given by 

X s = {v|r, € (P t(Ts))\ VT S G (33) 

where P^(Tfi) = {P2(75) + ar H Pt, &t h € R} is the space of quadratic polynomials enhanced 
by a cubic “bubble” function P 6 over 7*: for the pressure, we identify 

Y s = {. q\ Ti € Pi (T & ). Vr 6 € ( 39 ) 

We also introduce spaces of polynomial functions defined on the edges only, 

Qk = {t\s G P*(7), V 7 G £{T h )} f) 2. (40) 

where k identifies the order of the polynomial over the edge 7 . 

We now define two subdomain local spaces. First, for the velocity, the working and truth 
subdomain local spaces are given by 

U h {Th) = (P 1(T h ))\ (41) 

and 

v„(Th) = (vlr. € (pj(r,)) 2 , vr„ e n T „} f}x T „, (42) 

respectively, where we recall that 7 Zt h is the set of h- mesh elements that constitute 7//. 
We also define corresponding spaces which now include the incompressibility constraint, and 
define the spaces 

Z 5 (T h ) = {v € Us(T h ) I d T „(v,q) = 0 | V 9 € Ms(Th)}. (43) 
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where 


M h (T„) = P 1 (Th), 


(44) 


and 

M h {T„) = {q\T H e Pi(7fc), € -Rt h ], (45) 

are the local pressure spaces. 

Finally, we can define the associated global spaces with and without incompressibility 
constraint as 


V s = {v € A** | v| r „ € U S {T„)}, (46) 

and 

W s = {v € X m | v| Th € Zs(Tff), VT„ € T H j, (47) 

for 6 = H and 6 = h. In essence, the Us{T H ) and Z${T H ) are Neumann spaces over each 
T h , for which Vs and Hi are the corresponding global representations. Note that Z${Th) 
imposes the necessarv ’’global incompressibility constraint, on the velocity thanks to the 
discontinuous pressure approximation. 


4 Bound Procedure 

In this Section we present the hierarchical procedure to calculate the bounds. The three 
principal steps in this procedure are: calculation of the adjoint on the H- mesh; calculation 
of the hybrid flux on the H- mesh; and local Stokes solves to obtain the bounds. 


4.1 The H- Mesh Adjoint Calculation 

First, we solve the Stokes problem (l)-(2) on the working mesh. We look for (u h-Ph) € 
Xu x Yh such that 

a(w,Utf) - d(w,p H ) = ^(w), Vw € A'//, (48) 

— d(uH.q) = 0, V^€k//. (49) 

Second, we solve for the output adjoint. We look for (V , // 1 A^) € A h ^ such that 

a(v4, w ) - d(w . \%) = -(± £° l ’(w) + 2a(w,u w ) - f V (w)), (50) 

-d(^,q) = -(±f° P (<?)), V(w,g) € Xh x Yh- (51) 

Note that because we require H p (l) = 0 we can consider the zero-average space Y H , since 
solvability is ensured. Equivalently, (51) is in fact satisfied over the larger space in which T' 0 
in (39) is replaced by Regarding computational cost, we remark that (50)-(51) needs to 
be solved twice, once for each bound: ± refers to the pair of solutions required for the lower 
( + ) and upper ( — ) bounds. If a direct solver is used, only one LU factorization is required 
for (4S)-(51) — lh e Stokes operator is, in fact, the same, and only the right-hand side of the 
equations changes. 
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We now define a linear functional F ± (v, T, V) which represents the forcing term and the 
pressure term in the stress balance equations of (50); this functional is introduced mainly to 
simplify the notation. In particular, for any function J - in X" and V in >'°, we write 


F ± {v,?,V) = £ 4(v|r„ \T,V), Vv g X", (52) 

Th€Th 

where for all Tjj in 7}/, 

f t h ( w \th\?i'P) = ± £t„( v ) + a T H {f\T H - v) - fj w (v) - d T „{v,V). (53) 

We can now view the stress balance equations of (50) as 

2<i(w, u H ) = -F ± (w; A|). (54) 

We can also introduce a second fashion in which to re-express (54), that is 

B±(v , uh) = 0, Vv G Xfj. (55) 


Here, for any function Q in .V”, 

B-(y,Q)= E B f„M Vv G X", 

Th€ 7 ~h 

where for all Th in 7//, 

i? £f/ (w, £) = 2ar H (w,^|j H ) + F^(w;V»^, A ^), Vw € AV H . 


(56) 


(57! 


4.2 The F-Mesh Hybrid Flux Calculation 

The hybrid flux will appear in our Lagrangian as a Lagrange multiplier that enforces the 
subdomain continuity constraints. Recall that, for the Crouzeix-Raviart elements, we only 
need to impose continuity for the velocity components. Our procedure here is to calculate 
the hybrid flux by appealing to the broken space. To start, we have 

b(v,y ± ) — B ± {v, Utf), Vv G Vh, (58) 


that is, for all Th in Th-, 


E 

~it h €.£(Th) 


(T 1Th 

a T H 


[ V • y ± | £(7Tw) d 3 = B$ h {v, u h ), Vv G U h (T h )- 

Jit h 


(59) 


In [21] we present two different approaches to approximate the hybrid flux for quadratic 
elements based on earlier work for energy-norm estimators [17, 5]. These techniques are 
based on an initial approximation which is then corrected with a Pj term to ensure solvability. 
In addition, a higher-order (quadratic) term is included to improve accuracy; the latter is 
not required, but should give sharper bounds. We describe here only the most promising 
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approach of the two, in which the initial approximation is a P 0 approximation. For reasons 
of simplicity we will present the lower ( + ) bound hybrid flux calculation; the upper bound 
proceeds in a similar fashion. 

We first introduce y+ € (Qo) 2 , yj € (Qi) 2 and y+ € ( Q 2 ) 2 , which are different polyno- 
mial edge functions used in the hybrid flux approximation, y* = y 7 +y 7 +y 7 • The constant 
contribution to the hybrid flux, y+, does not present any new subtleties; it is obtained as in 
[20]. The linear correction, y+, is defined for each component by 

y+h = <0l(x) + (60) 

where al and are real coefficients to be determined, and (/^(x), 0^(x)) are linear edge 
functions constructed to be orthogonal to (yj H |-y 5 V 5 r f /I'/) [1T> 10]- The function <p>j- N is the 
restriction of the linear basis function associated with vertex n of T H to element T H . For the 
quadratic approximation of the hybrid flux we define 


where p 7 : 7 -4 R is the quadratic function uniquely defined by the conditions: 

J Pi^Th ds = 0- (6^) 

and . 

jf/>?d S = M, (63) 

where n in (62) refers to either of the two vertices of 7. We also introduce such that 
the function and associated with Th span P 2 (T//). 

Given y+, we calculate the coefficients 0 £ and thereafter following the procedure in 

[20]. To wit, we solve 


<Tf fvrjv* + y* + y*)^ 






(64) 


in which we exploit the fact that 

f A H (y*+yt+v*)^ = I = (^t+Oh 

J'y *''7 ** 


/h 


(65) 


since the quadratic function p 7 is orthogonal to the linear functions pj H . To calculate 
we then follow the procedure in [5], that is, we solve 


P t h 

T T h 


f &T H (y+ + y* + y+)ds — b^ h (p p Th , «»), 


(66) 


where y^~ and y^~ are now known. Details, in particular as regards solvability, may be found 
in [21]. 

To summarize, we first evaluate the non-conforming approximation, y+, to the hybrid 
flux, as in [20]. However, this approximation does not lead to solvability of the equilibrium 
equation (58). To ensure solvability, we solve N local systems, (64) one for each node of Th, 
to determine the constants of the linear contribution to y+, (60). Finally, we look for a 
quadratic contribution which leads to (66). Alternative approaches are described in [21]. 
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4.3 The h- Mesh Subdomain Neumann Problem 

Before we solve the subdomain problem, we must compute an adjoint 0* on the h- mesh. For 
sharp bounds, 0ft should be close to '0//- In addition, the adjoint 0^, must be continuous 
for all T h in T H to be a valid Lagrange multiplier. Finally, 0 ft must satisfy an equilibration 
equation if we are to obtain meaningful bounds — as we will discuss below. Therefore, for 

all Th in Th, we look for 0 ft £ U^(Th), such that 

<zr H (v,0f - 0#) - d,T H {v,pk) = 0, Vv £ Lt(Th), (67) 

-0r„(0*,<0 = -( ± C° P H {q)) Vq6Mk(TH), (68) 

where _ , 

U^(T h ) = {v|r H € U h (T H )\ v|^ r// = 0 //Ut w ^ v 7 ^ S{T h )}. (69) 

In effect, (69) is simply the affine manifold which imposes on 0± the H- mesh adjoint values 
0 * on the boundary of T H : it is important to note that, on dT H and dT h , the bubble function 
vanishes, so that the trace of j>g on 3T H is in U h {T H ) — the h - mesh subdomain space. In 
fact, (67)-(68) indicates that 0^ is an incompressible H l semi-norm projection of 0 W onto 
the fine-mesh; ph in (67) is a "dummy" variable (Lagrange multiplier) which is not used in 
the remainder of this work. Note that, if q = 1 in ( 68 ). then since 1 is in Mh(Th ), 

- c/t h ( 0 *, 1 ) = - / 0 h • n‘' TH ds = -d TH (V4- 1 ) = - ( ± © 1 )) (~°) 

from (51); recall that rT T » is the outward normal on 7 t h with respect to T H • The system 
(67)-(68) is thus solvable; note the issue of solvability does not arise in (67) because we do 
not have any Neumann problems — all boundaries are Dirichlet. 

For the local subdomain problem, we now look for £ i h(Tn)- for all Th in Th- such 
that 


2 «r H ( w - ) - dr„(w. ) 

= -( ± <?> 

) - 6r w (w) + a T „(0*. w) 


dr^w.Aft)- T. (?Th 

/ w-y ± \ Ehl 

. H> ds), Vw £ U h (T H ), 

(71 

-rr w ef(r w ) 

J 1T H 



UT H -q) 

= 0, Vq £ Mk{Tn)- 

(72 


To verify solvability, we take v = 1 in (59). The right-hand side of (71) then vanishes because 
1 £ Uh(Th) €: Uh{Th )■ Note that the construction of 0ft is also essential: the equilibrium 
equation (58) includes 0# , but in (71), ip h appears; however, since a{^ h , 1) = a(rp> H , 1) = 0, 
we are still able to ensure solvability. 

In a more compact notation, we can introduce two functions £ IF ft and 11^ £ Mk such 
that U±\t h = and n*|r H = *t h - £ Th, where W h satisfies 

2a(w,W^) - c?(w,fij) = -T ± (w;0*, A^) + 6(w,y ± ), Vw £ V h , 

-d(U±-q) — 0, Vg £ A//,. 
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(73) 

(74) 



We make several remarks. First, we note that the K systems for are completely decou- 
pled, leading to very efficient inversion compared to the original h - mesh problem of (6)-(7). 
This cost reduction is considerable, especially when decoupling the Stokes problem, which 
is a larger system with a larger bandwidth than an elliptic problem. An additional advan- 
tage is that each of these subdomain problems may be easily solved in parallel. Second, 
an additional constraint is introduced on the subdomain problems to impose the local in- 
compressibility constraint on U±. As we will see, this is not required by the bound theory, 
however we expect that it improves the accuracy of the bounds. By imposing the incom- 
pressibility constraint, which is more expensive, we look for the solutions to local Stokes 
problems instead of local Poisson problems. We have not yet investigated the latter. Third, 
note that we solve two (one for each bound) local Stokes problem to project the adjoint onto 
the h- mesh. The cost of this additional solve is small, especially if we use direct solvers, in 
which case onlv one LU decomposition is necessary for both the adjoint and the subsequent 
velocity calculations. 

Finally, we can now calculate the bounds as 

(sk)LB{H) = r} + , (75) 

and 

{s h )uB{H) = -T ] -, (76) 

where 

r/ ± = -a{Ut,Ui) ( 77 ) 

Note that the upper bound can be interpreted as the lower bound in which the output is 
multiplied by —1. 


5 Proof of Bounding Properties 

The proof of the bounding properties of r/* is based on classical quadratic duality theory 
[28, 19]. The key feature of our approach is the construction of a Lagrangian with a quadratic 
objective function and linear constraints such that, at stationarity, this Lagrangian evaluates 
to the output of interest. We first derive an “energy” equality that provides the stabilization 
in our Lagrangian. We take the test function in (4S)-(49) to be the solution (u^.p^) which 
yields 

a(\ih,u h ) - d{u h ,ph) = i N { u/,), ( 7S ) 

d{u hl q) = 0, V,? € Y k . (79) 

Note that (78)-(79) reduces to a quadratic form in U* — a(u h , u h ) - i N {u k ) = 0 — because 
the term d(uh,Ph ) is zero. For inhomogeneous Dirichlet boundary conditions, a boundary 
function would be introduced, as in [20], to directly obtain boundary conditions for the 
adjoint; the error formulation of [17, 25] can also be pursued. 
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By adding the output functional to the quadratic form (78) we obtain a function that 
reduces to s h = f° v (u^) + t° P {ph ) when (v = u A , q = Ph)- More precisely, we have 

± $h = min i±(° (v) ± £° P {q) + «(v, v) - £ N (v)) , (80) 

(v,?)e5 ' 7 

where 

a(^r , v) - d(n ,q) = £ N {fi ), V/i € A - /,, 

S={(v,q)e W h x n d(v, A) = 0, VA € Y h , • (81) 

6(v,t) = 0, Vt € Q 2 

The set of functions S is a singleton (v = U/,, q = p h ) equivalent to the solution of the 
Stokes equations (4S)-(49). Note we could replace If k with 1 h , which would yield decoupled 
Poisson rather than Stokes subproblems, as described in the previous section; we consider 
here the arguably more accurate choice 14 h - 

From a mathematical point of view, the solution to (80) is equivalent to finding the 
saddlepoint of a Lagrangian, : (v.q.fJ- ,A,t) € 44/, xfi x A/, xh x Q 2 

£ ± (v,^, pi , A, t) = ± C° l {v)±£° P {q) 

+a(v, v) — i N (v) + a{n , v) - d{fi ,q) - t N {n ) (82) 

— d(v. A) — 6(v, t). 

Inserting F ± (v; /x , A) from (53) and regrouping terms so that subsequent simplifications are 
more obvious, we can rewrite (82) as 

^(v,?,/*, A,t) = [-a(v,v) )] 

+ 2a(v,v) + F ± (v;/x , A) - 5(v,t)] 

+ [-d(/x ,q)±e° P {q)] . (83) 

Our first goal is to show that this Lagrangian evaluates to q ± of (77) for (v,?./x ,A,t) = 
•, Aj^y*), where represents any value in Y h . Proceeding, we obtain 

+ [2a(W*,Mf) + Af) - Hfay*) 

+ (S4) 

We immediately see that the first bracket of (84) equals r?*, where we recall that 

V ± = -a(U?,U?)-£*(j>t). (85) 

It remains to show that all the other terms in (84) vanish. We observe that the second bracket 
in (84) is, from (73), d(U±, f[±) which is zero thanks to (74). Finally, the last bracket in (84) 


15 



vanishes due to the construction of the adjoint, since we have imposed -d{ij) h , •)±f’° P (-) = 0 
in (68). 

We conclude that , 

rj* = C ± {ti£,-,^ k ,A£,y ± ). (86) 

It then follows from the classical quadratic linear duality theory that 

r/ ± < ±Sh, or r} + < s h < -if, (87) 

if . . ± 

£±(^,-,0*, A^,y ± ) = min £*^,-,0*, Aj,y ± ). (88) 

To demonstrate (88), we expand our Lagrangian (82) for v = U± + w, n = xjj k , A = 
A±, t = y ± , to obtain, 

£±(77^ + w.-,0^.A^, y*) = £ ± (^ ± .-,0ft .A^y*) 

+ 2a{w,U±) + F ± (w:0*. A*) - b(w.y ± ) + -d(0* , •) ± (° P {-) (89) 

+ a(w.w), Vw € Wh- 

We observe that all the terms linear in w (the first bracket) reduce to d( w, fi*) from (73), 

which vanishes since w € Wh C V h is incompressible. The terms -d(^ h .-) ± £° P (-) (the 
second bracket) also vanish thanks to the construction of the adjoint (67)-(6S). The re- 
maining term a(w,w) is positive semi-definite, which thus proves (SS) . Note that it is the 
energy equality which allows us to consider non-exact Lagrange multipliers and still provide 
non-infinite bounds. 

More precisely, to avoid meaningless bounds we need to verify that when minimizing 
our augmented Lagrangian we do not obtain oo. To this end, two main concerns must 
be addressed. First, solvability of (i3)-(i4) is essential. Without sol\abihty the teims on 
the mht-hand side could tend to infinity as the test function tends to infinity. Second, 
equilibration between — and ( Q ) is also essential, because these terms are not 
controlled by any quadratic stabilization. Because both of the above conditions are satisfied 
we are guaranteed non-infinite bounds. However, there is nothing in the presentation that 
proves that the bounds should be sharp. For the moment, we can suggest that, since A^ 
and y* are the saddlepoints of the H - mesh approximation to our Lagrangian. they should 
thus be close enough to the /i-mesh saddlepoint to yield good bounds. 

6 Optimal Stabilization Parameter 

In this section we present a procedure by which to improve the sharpness of the bounds. To 
this end, we introduce a positive real number, k, to scale our output s, and we look for the 
bounds to this scaled output. We then provide a procedure by which to calculate the optimal 
k, that is, the k that will yield the sharpest bounds; to be more precise, we maximize our 
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lower bound, and minimize our upper bound. A different but equivalent approach in which 
we scale the entire energy equality (78)-(79) is presented in R 1 in [19, 29], 

Our strategy to find the optimal k is to write all variables as linear functions in k, and 
then derive the bounds as a function of k. This procedure does not change the bounding 
theory and our bounds remain rigorous; indeed, our choices of Lagrange candidates are still 
valid even if the adjoint and the hybrid flux are decomposed into different contributions. 
The key is that these candidates must remain in the appropriate spaces, so some attention 
must be given to the boundary conditions. 

First, we decompose and A# as 

* ± o± * i± 

1>H = VH + K/t PH i 

A H - a h + kA h ■ 

where H € A’// satisfies 

= — (2g(w,u//) — £* v (w)), Vw 6 X H , (92) 

—d(ip H .q) = 0, V<? 6 Vff, (93) 

and ft^ € Xfj satisfies 

a{ftftw) - d(w,A)f) = -(±C° r ( w)), Vw€A' h , (94) 

-d(ftftq) = -(±f° P (7))- VqeY„. (95) 

Note that u//, the solution to (48), only appears on the right-hand side of the equation. In 
fact, in both equations the operator is identical, and we can take advantage of this fact for 

direct solvers. We now write tp^ as 

«± * o± *i± 

d’h = *Ph + K1 Ph 1 ( 96 ) 

which needs to satisfy, for each element Tjj on the H- mesh. 

a Tf/ (w, - rp H ) - dr w (w,p°) = 0, VweAT, (97) 

-d TH (i>T,q) = 0. V ? eL fl) (98) 

and 

ar w (w,^T - V'ff) - ^r H ( w .pt) = 0 Vw€.Vh, (99) 

-rfT,«T,9) = < 100 ) 

- o± - o± 

Based on similar arguments as in Section 5, the boundary conditions for rp h - tp H and 
ft* - ft* are homogeneous Dirichlet. These two set of equations are similar to (67)-(6S) 
in all respects; we force continuity of the adjoint across the subdomain boundaries, and we 
impose an incompressibility constraint in the interior of each subdomain. 


(90) 

(91) 
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We now present the k decomposition of the hybrid flux. First, we need to define two 
functions F 0± {v\T,V) and F x± ( v,J r ,'P). In particular, for any two functions T and V in 
X * and y , we define for all v £ X* 


F 0± (v;^,P) = £ \t h -,F,V), 

Th€'Th 

F l± (v f ^ P) = E ^(vlr.ri^n 


where 


*rJ(v;.F,P) = «r w (^v)-d Tw (v,n-^(v), 


Fj : t{v-,J r ,V) = ar H (^, v ) - d TH {v,V) ± £j h (v). 


Th 

Finally, we introduce 
As in Section 4.2. we solve, for all T H in T H . the following equations. 


y± = y 0± + K-y 1 *. 


( 101 ) 

( 102 ) 


(103) 

(104) 

(105) 


E 


'>T J 

(7 


* 0± 


ITlj 


T “ I v-y^U,, )ds = 2ar„(v,4)-f“(v;V-y.A?f), 

J~tT H 

Vv £ Uh(Th )• 


E °Tn f v-y^l^jds = -F^v^Ajfr, VvgW 

-kt h € 5(T w ) 


(106) 

(107) 


We can now solve the /;-mesh problems, 


and 


2a(w,u£ ± ) 

~d( w,nf) 

= -F 0± (w;^ ± .A° ± ) + 6(w.y 0± ). 

(108) 


-d{ UftS?) 

= 0. V(w,^) £ Xh x Yh. 

(109) 

2a(w,u] l ± ) 

-rf(w,ni*) 

- -^ 1± (w;^i ± ,Aj ± ) + 6(w,y 1± ), 

(110) 



= 0, V(w,<?) € X H x y' H . 

(111) 


We will not address solvability of (10S)-(109) and (llO)-(lll) as it follows our usual 
(see Section 4.3). 

Using the same derivation as in (77), the bounds can be expressed as 


^(k) = -- (afujj* 4- /cu^uf + KUfc*) + 

K \ 

= X( a{i ^x ± ) + t N ^l ± 


* 1± 
h 


K 


■2a(uf, a;*) - - Kafu^.u'*). 


proof 


( 112 ) 


(113) 
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Differentiating with respect to k , we find 


>,?(*) = i(a(uj*.u l ± ) + e‘'’{i>T))-a(i' h ± .u' k ± ), (114) 

>?«(«) = < i13 ) 


To optimize our bounds we require ri±(K m± ) = 0, which yields, 


K = 


\ 


■ 0± " 0± ' + i N wT) 


a(u k ,u h 


a(u 


i± -i±) 
h > u /i > 


( 116 ) 


To prove that k is a maximum, we proceed as follows: we first recall that 77 * is a lower bound 
to ±Sh. It follows that the terms in ^ must be positive so that our lower bound does not go 
to +oc as k decreases. These same terms also enter in the second derivative (and the radical) 
making the second derivative negative for all positive values of k (and the argument of the 
radical positive). The arguments are somewhat more transparent with the error formulation 
of [17]. 

We will now make some remarks concerning computational cost; we wish to show that 
we need onlv two subdomain solves rather than four to calculate the bounds for the optimal 
stabilization parameter k‘. It is obvious that the numerator is the same in both the upper 
and the lower bound calculations because it does not depend on the output functional. In 
addition, we can show from (94)-(95) and (99)-(100) that = ~i> h . Furthermore, we 

note that the right-hand side of (110) only differs by a sign when replacing rp h by -xjj h , 
which leads to u^ + = -u^ _ . Finally, because a(, ) is a symmetric positive-semidefinite form. 
a(u^ + , uj, + ) = ^(u^ - , ujf ), and thus the denominator of (116) is the same for both the upper 
and the lower bounds. From the above arguments, we obtain that k~ + = 

It follows that, in fact, we only need to perform two subdomain solves to compute our 
optimized bounds, just as in the non-optimized case. For clarity we summarize the relevant 
identities: u° h + = il 0 h ~, u\ + = -ujf, = rp° h ~ , and = -ip h ~. These identities also 
lead to an interesting property that the average of the bounds is not affected by k. The 
average of the bounds is given by 

\(V + ~T)-) = -^(a(u° + ,u° + ) + f v (^r) 

-o(u°~, u° h ~) - e N (rp° k )) 

-a(u£ + ,u][ + ) + a(ur^I _ ) 

+ b' v (^ _ ) 

(a(u) l + ,u)[ + ) -a(Ufc".ui - )) • (117) 

From the above identities we observe that the terms in ^ and ^ all vanish. After replacing 
the remaining and ^ h ~ by -u), + and , respectively, we obtain 

\{V + ~V~) = -2a(u \ + ,ul + )-C N {4>l + ), (US) 
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Figure 2: Velocity field solution for Th = T(h 0 , i)- 


thus proving the desired result. 


7 Numerical Results 

We present here results for the Stokes problem for a periodic domain (Figure 1) in which the 
flow is driven by a pressure gradient. The velocity field solution of this problem is shown in 
Figure 2 for the coarsest mesh, T ( h 0 ,i)- The triangulations investigated, T ( h 0 ,R), are uniform 
refinements of the coarsest mesh T (Ho , i) shown in Figure 3a. The //-meshes, T H , correspond 
to R = 1,2, 3, 4, 6, and the truth h - mesh corresponds to T h = T( Ho , 12 )', % is shown 

in Figure 3b. Note that, for all the refinement values of R considered, we satisfy X H C X h , 
as required by the theory. We shall denote the effective working-approximation element size 

associated with triangulation Th = T(h q< R) by H = l/R. 

Two outputs are investigated, as defined in (10): the flowrate, s (1) , and the lift force on 
the body, The test function * (€ H 1 ^) ® H^Vl)) used in the lift functional (13) is 

defined to be continuous and piecewise linear over the T H in T(h„,i) with 

X = 0, in Q \ f Y , 

x = 0, on r i5 i = {1,2, 3, 4}, 

x ■ e 2 = 1, on r 5 , 
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where fi' contains all the elements of T(h 0 .\) that have an edge on r 5 . Another choice of 
(incompressible) x ' s presented in [21], which yields almost identical results but at a higher 
computational cost. 

The objective here is to rigorously bound the output associated with To this end, 

different H - meshes can be exploited. It is obvious that the cost of the bound calculations 
increases as finer H - meshes are used, that is, as R increases for T(h 0 ,r )■ However, finer 
//-meshes also lead to sharper bounds because the adjoint and the hybrid flux are more 
accurately approximated. In fact, an adaptive procedure similar to [24] could be developed 
to efficiently produce a H - mesh and associated bound gap within a desired value. 

We can easily relate our hierarchical mesh procedure for calculation of the bounds to 
engineering design procedures based upon a hierarchy of numerical approximations. In fact, 
the first discretization, here the H- mesh, is a "working" coarse mesh approximation which 
is relatively inexpensive, but which generates solutions and associated outputs sg that are 
deemed sufficiently accurate for the purposes of "preliminary" analysis. The second dis- 
cretization. here the h- mesh, is a "truth" mesh which produces a solution and associated 
outputs s h for which \sh — s\ is assumed negligibly small. The /(-discretization serves to verify 
the prediction of the //-discretization, either prior to design, as in a validated— surrogates 
framework [32], during design, as in the trust-region optimization techniques [6], or after 
design, as final confirmation of the anticipated performance. Our bound procedure provides 
reliability of the truth mesh but at much lower cost. 

We plot in Figure 4a and 4b ($h)i:g/sh, )pref s hi (^OLb/^ an< ^ s H/$h as a function 
of (effective) H for, respectively, s (1 ) (flowrate), and s (2) (lift force). The average of the 
lower and upper bounds is denoted by (s^)* re . For the coarsest mesh, we observe that the 
upper bounds for both outputs are within ±15%. The accuracy of the lower bounds depends 
on the output considered. For the flowrate output. s (1 >, the lower bounds are within —5% 
and almost equal to s H \ in fact, in this case we have a "weak" compliance. (By compliance 
we refer to the property that the output calculated on the H - mesh is equal to the lower 
bound, {sh)LB = ?/ + = sg, which occurs when (i) the inhomogeneity of the weak form equals 
the output functional, (ii) the boundary conditions are homogeneous Dirichlet, and (iii) the 
operator of the problem considered is symmetric [21].) For s (2) , the lower bounds are within 
-20% of Sk calculated on T(h 0 . i)- ^e a ^ so observe that, for a refinement of two, both upper 
and lower bounds are well within ±10%. Recall that one of main advantages of the bounds is 
the certainty that s h does indeed lie within the calculated values. In practice, the “working” 
mesh should be constructed sufficiently accurate (considerably more so than our 7j// 0 ,i) used 
here). 

In Figure 5a and 5b we plot = log [(s/Jf-B ~ e LB = l°gl( 5 0LB — e pre = 
log \(s h )' - s*|, and eg = log \sg - s*| as a function of log H for s (1) and s (2) , respectively. 

For (sh)^ and sg appear to converge to Sh as 0{H 1 ' 5 ) as H — ► h. We would expect, for 
a smooth solution, that sg will converge at least as fast as 0{H 2 ), and no doubt faster. The 
corner singularities are most probably responsible for sg converging to s h only as 0(// 15 ). 
Note that from our “weak” compliance analysis in [21] the hybrid fluxes are zero, and we 
therefore rule out any error contribution from that calculation for the lower bound; as ex- 
pected, we obtain the same convergence rates for both {sh)l B and sg. Now, considering the 
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Figure 3 : (a) Coarsest working mesh T H = T [HoA) , and (b) truth mesh J H = T iHo , 12) • 
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e 4: Plots of (Sh)uB/ S h i ( 5 h)pre/^’ (■ s (i)lb/' s ^! an( ^ 5 
top) s (1) , the flow rate and (b, bottom) s (2) , the lifi 







convergence of ( $h)uB ( s< dll for the flowrate output), we note that we achieve only 0{H 1 ' 3 ) 
compared to 0(H 15 ) for ( s h )£ B ■ believe that this may be caused by the hybrid flux 

approximation — unfortunately preliminary work with a Pi initial approximation did not 
indicate any improvements [21]. Considering now s (2) , the quantities {s h ) m UEh ( s h )l B and 
( 5 0pre a N converge at the previous lower rate, 0(H 13 ), and s H converges at the same rate as 
for the flowrate, 0{H 1 ' 5 ). The same comments regarding the hybrid flux and the singularity 
can also be evoked for the lift output, s (2) , and no doubt explain the convergence rate results. 

The bounds presented here reflect the use of the scaling parameter k described in Section 6. 
For the flowrate output, n = 1 is optimal for all H (again due to the compliance result), 
while for the lift output k* tends to 0.0886 as R increases. Note that the choice of X does 
not influence significantly the accuracy and convergence of the bounds, as shown in [21]. 

We conclude with a few suggestions to improve the bounds for outputs of the Stokes 
problem. First, closer examination of the hybrid flux calculations is warranted; in particular, 
investigation of a P! initial approximation of the hybrid flux should improve the convergence 
rate of the bounds. Second, implementing the bounds technique within the stress formulation 
of the Stokes equations will allow for cleaner derivation of the lift linear functional. And 
finally, additional application to more relevant engineering problems will be presented in 
future papers. 
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the lift force. 
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